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Self-organization of moving objects in hydrodynamic environments has recently at¬ 
tracted considerable attention in connection to natural phenomena and living systems. 
However, the underlying physical mechanism is much less clear due to the intrin¬ 
sically nonequilibrium nature, compared with self-organization of thermal systems. 
Hydrodynamic interactions are believed to play a crucial role in such phenomena. To 
elucidate the fundamental physical nature of many-body hydrodynamic interactions 
at a finite Reynolds number, here we study a system of co-rotating hard disks in a 
two-dimensional viscous fiuid at zero temperature. Despite the absence of thermal 
noise, this system exhibits rich phase behaviours, including a fiuid state with diffusive 
dynamics, a cluster state, a hexatic state, a glassy state, a plastic crystal state and 
phase demixing.We reveal that these behaviours are induced by the off-axis and many- 
body nature of nonlinear hydrodynamic interactions and the finite time required for 
propagating the interactions by momentum diffusion. 


Self-organization is the autonomous organization of 
components into patterns or structures and the underly¬ 
ing mechanism can be grouped into two classes: thermal 
and athermal origin mm- Compared to the understand¬ 
ing of self-assembly of a thermal system, which can be 
explained in terms of the free energy, that of an athermal 
system is far behind due to the lack of a firm theoreti¬ 
cal background. Typical examples for the latter is ac¬ 
tive matter in hydrodynamic environments, which shows 
many interesting unconventional pattern formation with 
exotic dynamical features m [Mil]- Directional swim¬ 
mers undergoing translational motion, like fishes, tend 
to swim together while forming a cluster of particular 
shape mmm- The hexatic ordering of directionally self- 
propelling particles was also recently studied [12]. 

Rotation is another important type of motion. Self¬ 
organization of passive rotors was first observed in lab¬ 
oratory experiments by Grzybowski, Stone, and White- 
sides |4]. Such rotation can be induced by applying a 
torque to a particle by magnetic niniiii] and optical 
fields [15] . Self-organization of active rotors has also been 
studied intensively. For example, it was found that there 
are interesting stable bound states of spinning Volvox 
algae [8]. Furthermore, it was shown m that micro¬ 
organisms, like spermatozoa, self-organize into dynamic 
vortices and they form an array with local hexagonal or¬ 
der. This study indicates that large-scale coordination 
of cells can be regulated hydrodynamically, and chemi¬ 
cal signals are not required. It was also shown recently 
that self-assembly of rotors is a generic feature of aggre¬ 
gating swimmers m- There has also been a theoreti¬ 
cal prediction for intriguing self-organization of rotating 
molecular motors in membranes [18]. This and related 
problems have also been investigated by numerical sim¬ 
ulations [IHIISHSI]. 

As described above, there exist two types of rotors [25] : 
Rotors driven by external torques are called passive ro¬ 


tors, while those that are internally driven are called ac¬ 
tive rotors and many such examples can be seen in bio¬ 
logical systems. A realistic realization of a truly active 
system of self-rotors in biological systems may be one in 
which the particles are torque dipole with no resultant 
net torque on the system [H EH Eg. To elucidate the 
physics of self-organization of flow created by rotors at 
a finite Reynolds number, a system of hard disks each 
of which is rotated by an externally applied torque is 
an ideal model system. Ordering of this system may be 
regarded as a nonequilibrium counterpart of thermoday- 
namic ordering of hard disks. We note that , unlike a 
thermodynamic system, where a state is selected solely 
by free energy, dynamic factors such as hydrodynamic 
interactions also affects the selection of a state of an out- 
of-equilibrium system. 

For example, it has been now recognized that nonlin¬ 
ear hydrodynamic interactions play crucial roles in self¬ 
organization of rotating disks HEHIl]. In other words, 
the phenomena are beyond Stokes approximations, and 
may even have a link to self-organization of vortices [26] , 
which is observed at a high Reynolds number. Vortex 
crystals are very interesting such examples [271128]. The 
nonequilibrium, nonlinear, nonlocal, non-instantaneous 
nature of hydrodynamic interactions makes analytical 
approaches to this problem very difficult and thus numer¬ 
ical simulations are expected to play a crucial role. This 
problem also has a link to self-organization of a point 
vortex system [29l [30], but the finite size and solidity 
of disks lead to far more rich behaviour. Here we em¬ 
ploy a fluid particle dynamics (FPD) method [31], which 
we developed for studying hydrodynamic interactions be¬ 
tween colloidal particles (Methods). This method treats 
a solid colloid as an undeformable fluid particle inside 
which the viscosity is considerably higher than the sur¬ 
rounding liquid. This approximation makes us free from 
solid-fluid boundary conditions, which significantly sim- 
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plifies the computation. This method can quite naturally 
deal with many-body hydrodynamic interactions even at 
a high Reynolds number. 

In this Communication, we study self-organization of 
rotating hard disks in a two-dimensional (2D) incom¬ 
pressible liquid by using the FPD method. Like the Ising 
model for magnetic ordering or the hard sphere system 
for crystallization, this system may serve as a fundamen¬ 
tal model system for studying dynamical phase behaviour 
caused by hydrodynamic interactions between rotating 
particles. The situation is similar to the above-mentioned 
experimental and theoretical works mnnii]. We make 
our simulation in two dimensions (2D) to compare with 
the behaviour of its thermal counterpart, 2D hard disks, 
whose thermodynamic behaviour is reasonably under¬ 
stood [32]. Here we report surprisingly rich phase or¬ 
dering behaviours, such as aggregation, re-entrant order- 
disorder transition, glass transition, plastic crystal for¬ 
mation, and phase demixing in this class of strongly 
nonequilibrium systems. 


RESULTS 

Behaviours of a single and a pair of rotating disks 

Before discussing many-body interactions between 
disks rotating with an angular frequency D, first we de¬ 
scribe the behaviour of a dilute limit: Behaviours of a 
single rotating disk and a pair of rotating disks (Sup¬ 
plementary Fig. 1 and Supplementary Note 1). We 
characterize the rotation speed of a disk either by the 
angular frequency Q or the relevant Reynolds number 
Re = pa^Q./r]g {p: density; a: particle radius: pf. liq¬ 
uid viscosity). Figure and b show the 2D flow field 
around a rotating disk and the velocity distribution, re¬ 
spectively. The latter shows that the rotating velocity 
linearly increases with the distance from the centre of 
mass, r, inside a disk and decays as 1/r in its outside. 
This is the characteristic of the so-called Rankin vortex, 
i.e., a forced vortex in the central core surrounded by a 
free vortex. We note that the Rankin vortex is known 
to mimic the 2D flow field of tornado and hurricane [33| . 
For two co-rotating ‘point’ vortexes, the analytical solu¬ 
tion is known and the two particles rotate around their 
center of mass towards the same direction as the rotating 
direction [HEn]. In this case, the radius of rotation is 
the half of the initial interparticle distance. For particles 
with a finite size, on the other hand, the direction of ro¬ 
tation and the rotation center are the same as the case of 
point particles, but the radius of the rotation decreases 
with an increase in the rotation speed of the particles 
U, or Re. This problem has a similarity to viscous in¬ 
teractions of co-rotating vortices 1211 IS], but the crucial 
difference arises from the fact that the vortex core is an 
undeformable solid and not a liquid in our case. There 
are a few studies on spinning particles in three dimen¬ 
sions [23 ESI; however, we note that there is an essential 


difference between 2D and 3D problems (Supplementary 
Figure 2 and Supplenetary Note 2). We confirm that if 
we fix the centres of mass of two rotating particles on a 
fixed line, they always repel with each other by the Mag¬ 
nus force HlHEll. The attraction between particles is 
thus due to interparticle hydrodynamic interactions, as 
schematically explained in Fig. [^. The flow field gener¬ 
ated by a rotating disk makes the other disk follow it, and 
vice versa. Thus, two particles try to follow each other 
while rotating around the centre of mass of the pair. In 
3D, on the other hand, such hydrodynamic interactions 
are weaken due to the presence of the escape dimension, 
thus the Magnus force wins over the hydrodynamic at¬ 
tractive force, and particles repel with each other. In 2D, 
this hydrodynamic attraction overwhelms the repulsion 
due to the Magnus force, which leads to the rotating pair 
of particles (Fig. Hi), whose interparticle distance mono- 
tonically decreases with an increase in D, or Re (Fig. 
e), as long as the area fraction of disks, is sufficiently 
small. However, it should be noted that this situation is 
realized in a periodic boundary condition. In relation to 
the above dimensionality effect on hydrodynamic interac¬ 
tions between rotors, it is worth noting that the hexago¬ 
nal ordering observed by Grzybowski et al. [DiiiiaiTii at 
finite Reynolds numbers is due to the effective hydrody¬ 
namic repulsion between the disks floating at an interface 
in three dimensions (z.e., 2.5 dimensions). 


Structural ordering due to many-body 
hydrodynamic interactions 

Now we consider the dynamical behaviour of a system 
of many disks rotating in the counter-clockwise direction 
with O. The hydrodynamic attraction between rotating 
disks leads to the formation of rotating clusters for low $ 
(Fig. [^). The higher rotation speed of individual disks 
leads to the formation of a more compact rotating clus¬ 
ter. This tendency is basically the same as that for a pair 
of rotors (Fig. El)- For this regime, only one cluster is 
formed in the simulation box and the whole cluster ro¬ 
tates in the counter-clockwise direction, as shown in Fig. 
[^. A cluster always tends to have a circular shape, but 
it does not have any particular internal structural order, 
partly because imperfect matching between the size and 
the number of particles leads to structural fluctuations: 
the internal structure is basically controlled by the num¬ 
ber of disks in it and Re, which determine the cluster 
size, but fluctuating with time. With a further increase 
in Re, however, this cluster state becomes unstable, since 
the repulsive Magnus force of nonlinear origin eventually 
wins over the hydrodynamic attraction for high Re (see 
below). Above a critical ^ (^c ^ 0.05), on the other 
hand, a system exhibits a re-entrant transition between 
states as a function of D (Fig. |^): for low D a system is 
in a disordered liquid state with large fluctuations, but 
with an increase in Q it enters into a rather stable hexatic 
phase where rotating particles are localized on a hexago- 
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FIG. 1. Behaviours of an isolated rotating disk and co¬ 
rotating disks in a dilute suspension, a, 2D flow fields 
around an isolated rotating disk, b, The velocity profile as 
a function of r/a. v/[au) increases linearly with r inside 
the disk and then decays as 1/r in its outside, which is the 
characteristic of the so-called Rankin voltex. c, Schematic 
picture of a pair of two disks rotating around the centre of 
mass. The two disks are rotating counter-clockwise with the 
same speed, d, Trajectory of coupled rotating particles (the 
system size=256x256). The initial separation of the particles 
(arrowed) was 30 and each particle rotates counter-clockwise 
with Re = 0.76. e, Temporal change in the interparticle sepa¬ 
ration. We note that for all the cases including Re = 0.7, the 
interparticle distance eventually reaches a final steady-state 
value of the separation (see Fig. [^), which monotonically 
decreases with Re. f, A cluster of disks formed at $ = 0.04 at 
Re = 5.94. We confirm that irrespective of the value of Re, a 
cluster is always formed in the range of Re studied. 


nal lattice. The border between the cluster to the hexatic 
state is rather sharp as a function of 

The transition can be characterized by the nature of 
interparticle interactions. Here we analyse a point pat¬ 
tern to extract the nature of interparticle interactions. 
The point pattern analysis is very useful to determine 
the overall interparticle interactions [37]. Here we use 
what we call N function, which is the number of con¬ 
nected regions as a function of the radius of circle whose 
centre is located at the centre of each disk. The num¬ 
ber of connected regions decreases monotonically, with 
an increase in the circle radius, from the total number of 



FIG. 2. Re-entrant state transitions observed at a 
rather high 4>. a, Snapshots of particle configurations to¬ 
gether with the velocity fields. The system size is 256x256, 
the number of particles is 80, and 4> ~ 0.157. We can see se¬ 
quential transitions from disordered, hexatic ordered, to dis¬ 
ordered states with an increase in Re. b, Upper panel: The 
direction of the hydrodynamic force 0 measured from the in¬ 
terparticle axis. Lower panels: Patterns formed by Brownian 
simulations for the two off-axis forces: Left is a liquid state for 
0 = 89°, whereas the right panel is a hexatic state for 0 = 85°. 
c, Re-dependence of the hexatic order parameter Te, which 
clearly shows the re-entrant behaviour for various 4>’s (ex¬ 
pressed in %). d, The degree of fluctuations of Te as a func¬ 
tion of Re for various 4>’s (expressed in %). We can see that 
as in a thermodynamic hexatic ordering transition, the or¬ 
der parameter exhibits large amplitude fluctuations near the 
transition points, e, Re-dependence of the decay of the spa¬ 
tial correlation function of the hexatic order parameter g6{r) 
normalized by the radial distribution function g(r) around 
the transition at low Re for 4> = 0.157. In the hexatic state, 
gQ{r)/g(r) decays with a power law with the exponent of - 
1/4, as it should be [32]. In the disordered states it decays 
almost exponentially, f, The same as e but for around the 
transition at high Re. In the hexatic state, gQ{r)/g(r) again 
decays with a power law with the exponent of -1/4. In the 
disordered states it decays almost exponentially. 


disks Nq to one. Here we use N for the one normalized 
by Nq. By comparing N for a reference system made of 
randomly distributed disks, i.e., Poisson pattern, we can 
judge whether the interaction is repulsive or attractive. 
For a system of particles with repulsive interactions, par¬ 
ticles form a rather regular pattern and N decays slower 
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than that for the corresponding Poisson pattern. For 
a system of particles with attractive interactions, on the 
other hand, particles form a cluster pattern and N decays 
faster than that for the corresponding Poisson pattern. 
In Fig. [^, we show the N function for 4> = 0.04. We can 
clearly see that for the cluster-forming case at i?e = 0.72 
the interaction is attractive whereas for the disordered 
state at Re = 5.97 the interaction is repulsive. In Fig. 
[^, we show the N function for 4> = 0.15. We can see 
that for all the value of Re the interaction is basically 
repulsive, but its strength is maximum at Re = 5.9, for 
which the hexatic order is formed. 




FIG. 3. The analysis of point patterns, a, The decrease 
of as a function of the particle radius for 4> = 0.04, 

where p is the particle number density and is the aver¬ 

age interparticle distance, b, The decrease of A/" as a function 
of the particle radius ajp^^"^ for 4> = 0.15. 

The transition can be seen even more clearly by the 
degree of localization of the flow field: for a hexatic 
state each particle has its own localized rotational flow 
field, whereas for a cluster state a single vortex is always 
formed and thus the flow field is strongly delocalized. 
However, the precise nature of the cluster-hexatic phase 
transition, such as whether the transition is continuous 
or discontinuous and whether a disorder state always ex¬ 
ists between the two states or there exists a critical point 
between the cluster and the hexatic state, is not clear at 
this moment. To access this problem, we need to survey 
the border region of a bigger system size with a high res¬ 
olution of 4> and Re. Although this is a very interesting 
problem, we leave this for future investigation. 

We also find that the hexatic state is eventually desta¬ 
bilized and melts by a further increase in ft and a system 
becomes disordered again. We note that for all these 
states interparticle interactions are basically repulsive, 
unlike the case of low 4>. Although the torque exerted 
to each disk is exactly the same, the rotation speed of 
a disk can in principle depend on the particle configura¬ 
tion around it. Here we show in Fig. |^the normalized 
variance of the angular frequency ft of rotating disks as 
a function of the averaged ft. For ordered states, the 
variance, i.e., the fluctuations of rotation speed, becomes 
very small, indicating all the disks rotating with almost 
the same frequency. For disordered states, on the other 
hand, the variance is large, reflecting the large fluctu¬ 
ations of particle environment. This result indicates a 
strong negative correlation between the degree of fluctu¬ 



FIG. 4. The normalized variance of the angular fre¬ 
quency Q of rotating disks as a function of the aver¬ 
aged Q for 4> = 0.15. We can see the negative correlation 
between the degree of the hexatic order and the variance. 

ations of rotational speed of particles and the degree of 
order 

In Fig. [^, we show the f^-dependence of the hexatic or¬ 
der parameter 4^6, which clearly indicates the re-entrant 
nature of the state transitions. The transitions can also 
be characterized by the magnitude of fluctuations of 4^6, 
or the susceptibility (see Fig. |^), which is also observed 
in a thermodynamic hexatic ordering in 2D disks. In the 
hexatic ordered state, we confirm the power law decay 
of the spatial correlation of the hexatic order which is 
specific to the hexatic phase (see Fig. and f). We do 
not see any indication of the positional order since the 
radial distribution function decays almost exponentially 
(see below). 



FIG. 5. Hexatic ordering in a large system. Here we 
show results of large-size simulation (lattice size=2048^ and 
the number of particles=5120). The area fraction 4> = 0.15 
and Re = 5.9. a, The structure of a hexatic state shown with 
three colours (yellow for particles with six neighbours; red for 
particles with more than 7 neighbours; blue for particles with 
less than 5 neighbours), b, The same as a but with orienta¬ 
tions of hexatic order. See the colour bar on the meaning of 
colour. 

Here we show a hexatic phase observed in a large 
system (lattice size=2048^ and the number of parti- 
cles=5120) at 4> = 0.15 and Re = 5.9. We can see grain 
boundaries between hexatic order with different orienta¬ 
tions in Fig. IT and b. Figure and b show the decay 
of the correlation function of the hexatic order normal- 
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ized by the radial distribution function g{r)^ 
and that of g{r) — 1, respectively. We can see gQ{r)/g{r) 
decays algebraically in a short distance r < 400, but de¬ 
cays faster for long distance. This is because the size of 
mono-domain regions is finite (see Fig. §1 and b). On 
the other hand, we can see that g{r) decays faster than 
an algebraic decay even for r < 400, where gQ{r)/g{r) de¬ 
cays algebraically. This clearly indicates the absence of 
quasi-long-range translational order in the ordered phase. 
Thus we conclude that the ordered phase is the hexatic 
phase. The appearance of the transitions between dy¬ 
namical states as a function of 4> and 91 in athermal sys¬ 
tems is quite striking, which is reminiscent of phase tran¬ 
sitions in a thermal system. We stress that the interpar¬ 
ticle interaction in our system is of purely hydrodynamic 
origin. 



FIG. 6. Spatial correlation of hexatic and positional 
order, a, The spatial decay of the correlation function of 
the hexatic order gQ{r) in Fig. normalized by the radial 
distribution function g{r). We can see it decays algebraically 
in a short distance, but decays faster for long distance because 
of the finiteness of the domain size. The dashed line has a 
slope of —1/4. d, The spatial decay of the radial distribution 
function g(r) for the pattern shown in Fig. It decays more 
quickly than the decay of ge (r) even for rather short distance 
(r < 400). The dashed line has a slope of —1/3. 


Here we consider a quite interesting feature of hydro- 
dynamic interparticle interactions, which are absent in a 
thermal system. Binary interaction potentials that lead 
to phase ordering in thermal systems always act along the 
interparticle axis, i.e., along the line connecting the cen¬ 
tres of mass of two interacting particles. However, this 
is not the case for our athermal system. The nonlinear 
Magnus force acts along the line connecting two particles, 
whereas the Stokes force acts according to the flow di¬ 
rection. Furthermore, linear hydrodynamic interactions 
are tensorial and act not along the interparticle direc¬ 
tion. Thus, the total hydrodynamic force does not act 
along the interparticle axis. The direction of the hydro- 
dynamic force as a function of Re is shown in the upper 
panel of Fig. with a schematic explanation. With an 
increase in i?e, both the Magnus and the hydrodynamic 
force increase. However, the nonlinear Magnus force in¬ 
creases more rapidly than the linear hydrodynamic force. 
Thus, the direction of the total force, which is repulsive, 
approaches the interparticle axis. A strong enough re¬ 
pulsive force acting on particles not so far from the inter¬ 


particle direction leads to the formation of hexatic order. 
To verify this scenario, we have performed Brownian dy¬ 
namics simulation (without hydrodynamic interactions) 
by changing the angle 0 between the direction of an arti¬ 
ficial repulsive force and the interparticle direction. We 
find that when 0 is continuously decreased, a system in¬ 
deed forms hexatic order below a critical value of 0 (see 
the lower panels of Fig. [^). 

Next we discuss the nature of the translational motion 
of rotating disks in the disordered states. To see this, we 
calculate the mean-square displacement of disks (Ar^(t)) 
(see Fig. Interestingly, in the two types of disordered 
states particle motion is apparently diffusional: In the 
long-time limit we observe the relation (Ar^(t)) ^ 
where I^eff is the effective diffusion constant and t is the 
time duration. It should be noted that in our system 
there is no thermal noise and an isolated rotating disk 
exhibits no motion. For a pair of rotating particles, we 
also observe their trajectory is very stable and there is 
no fluctuation (see Fig. W)- Thus, the apparently diffu¬ 
sional motion should be the consequence of self-generated 
force noise due to many-body hydrodynamic interactions 
of both linear and nonlinear origin. We can also see no 
diffusion, or non-ergodic behaviour, for the hexatic state. 
The random nature of fluctuations may be related to the 
long-range nature of hydrodynamic interactions, which 
makes a number of the surrounding particles affecting a 
particle large enough to provide strong stochastic spatio- 
temporal fluctuations. 

For high i?e, the particle diffusion starts to become 
comparable or faster than momentum diffusion: Deff ^ 
z^, where u is the momentum diffusion constant or the 
kinematic viscosity v = gjp (see Fig. and b). This 
implies that hydrodynamic interactions induced by the 
rotation of a particle cannot fully propagate to its neigh¬ 
bouring particles for high Re. This weakens the repul¬ 
sive interactions of particles and eventually leads to the 
melting of the hexatic state. Thus the re-entrant hexatic 
ordering as a function of Re may be explained as follows: 
The increase of the Magnus force of nonlinear origin with 
an increase in Re makes the direction of the interparticle 
force more aligned along the interparticle direction and 
also increases the strength of the repulsive force. This 
leads to stabilization of the hexatic order, as explained 
above. The ordered state is stable until the repulsive in¬ 
teraction is weakened by the intrinsically kinetic nature 
of hydrodynamic interactions: Unlike ordinary interpar¬ 
ticle interactions which propagate with the speed of light, 
hydrodynamic interactions propagate much slower in a 
diffusive manner with the momentum diffusion constant 
v. This kinetic weakening of the repulsive force eventu¬ 
ally destabilizes the hexatic state and leads to the melting 
into the disordered chaotic state. 

Here we summarize what we observed in our system 
and show the state diagram as a function of 4> and Re 
(see Fig. ^). We can see the three states, i.e., disor¬ 
dered fluid, cluster, and ordered hexatic state. It is quite 
striking that particles interacting by hydrodynamic in- 
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FIG. 7. Dynamical behaviour and state diagram of rotating disks, a, The mean-square displacement vs. 

time t. Despite the absence of thermal noise, particles undergo Brownian-like diffusive motion = Dgfft with an effective 

diffusion constant Deff in the long-time limit. The black dashed line is the momentum diffusion constant, or the kinematic 
viscosity z/(= T]i/p)- b, i^e-dependence of the effective diffusion constant Deff. In the hexatic ordered state, Deff is very low 
and there is no diffusive behaviour, as it should be (see a). At high Re, Dgff exceeds u (indicated by the blue horizontal line), 
which means that fluctuations of particles cannot be suppressed by hydrodynamic interactions at high Re (see text), c, State 
diagram on the Re-^ plane. At low $ and low Re, the system forms a cluster due to hydrodynamic attractive interactions. 
At high on the other hand, the system exhibits re-entrant transition between a disordered chaotic liquid state and ordered 
hexatic state. The latter state is basically stabilized by the repulsive interaction due to the Magnus effect. 


teractions alone exhibit such rich phase behaviours. 


Other interesting states formed by rotating disks 

Finally, we show other interesting states formed by ro¬ 
tating particles. The introduction of size polydispersity 
to hard spheres is known to lead to the formation of a 
glass state for a thermal system [38]. Motivated by this, 
we introduce the polydispersity in the rotating speed of 
particles, whose variance is 6, and indeed find a non- 
ergodic glassy state of the rotating particles (Fig. If) 
for high enough 5 (Fig. [^): liquid-glass transition in an 
athermal system. Reflecting the nearly continuous na¬ 
ture of the liquid-to-hexatic transition [32l [39] , there is 
no sharp transition from a hexatic to a disordered glassy 
state. Although disks are rotating around their centres of 
mass, their positions are frozen in a disordered configura¬ 
tion and thus the system can be regarded as a nonergodic 
glassy state. Here we show how the distribution of the 
rotation speed of particles leads to the loss of hexatic 


order and results in the formation of a glassy state. As 
shown in Fig. If, the larger deviation from the aver¬ 
aged Re leads to the larger deviation from the average 
number of nearest neighbours (=6). With an increase 
in the variance of the distribution 6, more defects are 
produced and the hexatic order is eventually lost above 
6 > 0.2. We can also see that the the interparticle dis¬ 
tance is larger for a particle rotating with a faster speed 
because of stronger hydrodynamic repulsive force, i.e., 
Magnus force (Fig. [^. Both of these disorder effects 
are responsible for the formation of a glassy state: The 
number of the nearest neighbours and the average dis¬ 
tance to the neighbours are both strongly correlated to 
the rotational speed of particles, which explains a wide 
enough distribution of the rotation speed results in the 
formation of a non-ergodic amorphous state instead of a 
hexatic ordered state. This glassy state is non-ergodic 
if we consider the particle configuration, yet maintains 
strong flow fields, which makes this state very unique. 
Reflecting the kinetic origin of interparticle interactions, 
the introduction of disorder in the dynamic quantity, Vl, 
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FIG. 8. Other interesting states formed in a system made of rotating particles. a, A glassy non-ergodic states 
formed in a system of rotating disks where the torque F acting on particles has a Gaussian distribution whose normalized 
variance is (5 = (AF^)^/^/F. Here $ = 0.157, the average Re {Re) = 4.22, and S = 0.2. The colour of particles is green when 
the number of nearest neighbour particles NN is 6. For NN > 6, the colour is red and for NN < 6 blue, b, (5-dependence of the 
hexatic order parameter Te. Here = 0.157 and Re = 7.76. With an increase in 6, the hexatic order monotonically decreases 
and the system eventually enters into a nonergodic glassy state, c, A plastic phase formed in a system made of rotating 
dumbbells (counter-clockwise). The torque applied is the same for all particles. The rotation direction is counter-clockwise. 
The system exhibits hexatic order at a certain range of Re, but without any orientational order of the axes of dumbbells, d, 
fe-dependence of Te for rotating dumbbells. We note that the transition is rather broad. The schematic image represents 
the distribution of volticity around a single rotating dumbbell (red: counter-clockwies; blueiclockwise). e, Phase separation of 
particles rotating clockwise (Q) and counter-clockwise (—(f^ = 0.155 (or. Re = 6.35)). 4> = 0.157 and the number fraction 
of particles rotating clockwise is 0.5. 
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FIG. 9. Correlation between local structure and rota¬ 
tion speed (or, Re), of each particle in a glassy state. 

a. Dependence of the number of nearest neighbours on Re for 
individual particles. The system size is 256^, (j) = 0.15, and 
the average value of Re, Re, is 4.2. We can see that the larger 
deviation from Re leads to the larger deviation from the av¬ 
erage number of nearest neighbours (=6). b. Dependence of 
the distance to nearest neighbour particles on Re for individ¬ 
ual particles. The conditions are the same as the above. We 
can see that the the interparticle distance is larger for a par¬ 
ticle rotating with a faster speed (or, larger Re) because of 
stronger hydrodynamic repulsive force, i.e., Magnus force. 


is essential for avoiding the ordering, which is an inter¬ 
esting point unique to purely kinetic athermal systems. 

We also find a plastic crystal-like state in a rotating 
particle pair (dumbbell) system (Fig. |^) above a criti¬ 
cal Re (Fig. ^). This state is characterized by hexatic 


ordering of the centre of mass of disk pairs (dumbbells) 
without any orientational order in the axis directions of 
dumbbells. So we find for this type of athermal systems 
made of rotating particles almost all states seen in its 
thermodynamic counterparts, including a liquid, a hex¬ 
atic phase, a glass, and a plastic crystal. The only miss¬ 
ing state is a liquid-crystalline state, but the lack of this 
state is natural consequence of the fact that the system 
is composed of rotating elements. We also note that the 
modification of the rotational direction leads to complex 
behaviours; for example, we can introduce phase demix¬ 
ing of particles rotating oppositely (Fig. [^), which may 
be used to separate different types of passive and active 
rotors. In relation to this, it is worth noting that phase 
separation between particles rotating clockwise and anti¬ 
clockwise was recently observed even without hydrody¬ 
namic interactions m- It is interesting that such phase 
separation is observed in both systems with and without 
hydrodynamic interactions. It is also worth noting that 
the presence of regular arrays of vortices, e.g., the trian¬ 
gle state, are predicted for an active polar film [41]. So 
far we have not seen such a regular state, but the similar¬ 
ity in the physics between the two systems implies that 
it might exist in a certain parameter range. This is an 
interesting problem for future study. 













DISCUSSION 

It is remarkable that our athermal system where par¬ 
ticles interact only via hydrodynamic interactions ex¬ 
hibits such rich phase (or more strictly, state) behaviour 
and reproduces almost all the physical states observed in 
its thermal counterpart. We hope that these phase be¬ 
haviours will be observed experimentally. For this pur¬ 
pose, a quasi-2D version of the experimental setup used 
in [DiiiniiTii may be suitable, since the 2D nature of hy¬ 
drodynamic interactions is important. Here we have fo¬ 
cused on the bulk behaviour of rotating disks to study the 
fundamental nature of dynamical phase behaviour from 
a viewpoint of nonequilibrium statistical physics. The 
effects of confinement on such a system are also quite in¬ 
teresting not only from a fundamental viewpoint, but also 
from a viewpoint of applications to microfluidics [42j |43] . 
Such effects on rotors in a fluid have recently been numer¬ 
ically studied by Gotze and Gompper [22l [23] . It is quite 
interesting to study how such spatial confinements affect 
all the dynamical phases we reported and the transition 
between them. 

Here it is worth stressing that even non-ergodic states 
of particles such as hexatic, glassy, and plastic crystal 
states have strong hydrodynamic flow fields, which makes 
these states quite distinct from their thermal counter¬ 
parts. The interesting and unique feature of hydrody¬ 
namic self-organization is that structural ordering is the 
consequence of self-organization of flow dissipating en¬ 
ergy and thus even non-ergodic states are maintained by 
dynamical flow. In nature, there are many dynamical 
systems in which crucial interactions between elements 
are of purely hydrodynamic origin. The coexistence of 
linear and nonlinear hydrodynamic interactions, the re¬ 
sulting unconventional off-axis force, the finite propaga¬ 
tion speed of the interactions, and the significance of hy¬ 
drodynamic degrees of freedom even for non-ergodic (ap¬ 
parently static) states lead to rich and non-trivial self¬ 
organization. We hope that our study sheds new light 
on hydrodynamic self-organization and stimulate further 
study on this intriguing problem. 


METHODS 

Simulation method. Treating hydrodynamic inter¬ 
actions between colloids is difficult even for a thermal 
system. There are several methods such as Stokesian 
dynamics, lattice-Boltzmann, smooth-particle, and fluid- 
particle-dynamics (FPD) methods. Here we employ our 
FPD method [HI [44], which has an advantage in its 
theoretical transparency and its applicability to a high 
Reynolds number {Re) regime. We can access a high Re 


regime rather easily particularly because a torque applied 
externally makes the flow field inside a fluid disk almost 
exactly that for a solid disk even at high Re. 

Here we briefly explain the FPD method [H] and the 
physical concept behind it. A particle whose centre of 
mass is located at is represented by a smooth viscosity 
change as 7?(r) = rje + {tJc - rit)4>i{r), where % is 
the liquid viscosity and r]c is the viscosity inside a colloid 
particle. The summation is taken over all N particles. <fi 
represents particle i as (j)i{r) = {tanh[(a — 1“^ — '^i|)/^] + 
l}/2, where a is the particle radius and f is the interface 
thickness. Then the equation of motion to be solved is 

p(^ + 'y-v)-w = /t/ + /T-v-n (1) 

with n = pi — ? 7 (Vi;i ^ Vv), where I is the unit 
tensor. Here p is the mass density, and we assume 
that the density of the liquid is the same as that of 
particles. v{r) is the velocity field, and the pressure 
p is determined to satisfy the incompressibility condi¬ 
tion V • i; = 0. Here is the force density due 

to the interparticle interaction determined as fu{r) = 

where U{r) is the 

interparticle potential, Vij = Vi — r^, and A = Ai = 
f dr(j)i{r) is the area of each particle, /^(r’) is the force 
density due to the torque: /^^(r) = a\r\(j){r)eo^ where 
ee is the unit angular vector in the counter-clockwise di¬ 
rection. a is the strength of the torque and a > 0 leads 
to the counter-clockwise rotation of a particle. 

In our FPD method the particle rigidity is approx¬ 
imately expressed by introducing the smooth viscosity 
profile, ri{r). The approximation is better for a larger vis¬ 
cosity ratio ric/ri£ and a smaller f/a. By multiplying both 
sides of Eq. 0 by (j)i{r) and then performing its spa¬ 
tial integration, we can straightforwardly obtain an ap¬ 
proximate equation of motion of particle i\ MidVi/dt = 
Fi + Ki^ where Mi = pA = M and Vi = f drv(pi/A are 
the mass and the average velocity of particle i, respec¬ 
tively. On the right hand side, Fi = f dr(j)i{fjj + /^) 
is the force arising from the interparticle interaction and 

the torque, and Ki = — f dr(j)iV • H = — / dSihi • H the 
force exerted by the fluid. Here we use the following ap- 

■f-)- "f-)- 

proximate relation j drV(l)i - Q = — riidSi • Q for an 

•f-)- 

arbitrary tensor Q(r), where Si is the surface of particle 
i and hi is the unit outward normal vector to Si . In prac¬ 
tical numerical calculations, the on-lattice velocity field, 
v{r^ t + At), is evaluated from the physical quantities at 
time t by Eq. 0 . Then we move particle i off-lattice as 
a rigid body by Vi {t + At) = Vi (t) + AtU i (t + At), where 
At is the time increment of the numerical integration. 

In our simulation, the units of length I and time r 
are related as r = j {rnjp)., which sets both the scaled 
density and viscosity of the fluid region to unity. This 
r is a time required for the fluid momentum to diffuse 
over a lattice size i. The units of stress and energy are 
G = p{i/T)‘^ and e = respectively. Eurthermore, we 
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set Ti^jrig = 50, At = 0.003, i = 0.5^ = 0.5, and a = 6.4. 
We confirmed that this choice yields reliable results by 
comparing them wit the analytical solution for a single 
particle rotation (see, e.g.. Fig. and b). The simula¬ 
tion box used was typically = 256^. To avoid cumber¬ 
some expressions, we will use the same characters for the 
scaled variables below. We solve the equation of motion 
[Eq. Q] by the Marker-and-Cell (MAC) method with a 
staggered lattice under the periodic boundary condition. 
Interparticle potentials. To mimic rotat¬ 

ing hard disks, we employ the Weeks-Chandler- 
Andersen (WCA) repulsive potential [45]: Ujk{r) = 
4e - (cTjfe/r)® + 1/4} for r < 2^ajk, oth- 

erwise Ujk{r) = 0, where e gives the energy scale, 
ajk = {<Jj + cr/c)/2 and aj represents the size of particle 
T 

Characterization of structures. The 2D radial dis¬ 
tribution function g{r) was calculated as 


which is the ratio of the ensemble average of the number 
density of particles existing in the region r ^ r + Ar 
to the average number density p = Njl?. Here N is 
the number of particles in the simulation box, whose side 
length is L, and Ar is the increment of r. 

Similarly, the spatial correlation of Tg is calculated as 




The spatial correlation of the bond-orientational order 
can then be characterized by gQ^{r)lg{r). 
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SUPPLEMENTARY FIGURES 



Supplementary Figure 1. Relation between the 
exerted torque Teff and the angular velocity Q of 
a disk. 



Supplementary Figure 2. Switching of the effec¬ 
tive hydrodynamic interaction between 3D and 

2D. Decreasing the size of the third ^ dimension of the 
quasi-2D simulation box (128^ in x and y directions) from 
32 to 16 leads to the switching of the interaction from re¬ 
pulsive to attractive, or the switching of the relaxation 


mode of vortex from cascade to anti-cascade mode. Ac¬ 
cordingly, the two rotating disks can have a stable rotat¬ 
ing trajectory with a constant radius for the latter. 


SUPPLEMENTARY NOTES 


Supplementary Note 1: The behaviour of a single 
rotating disc 

As we show in the main text (Fig. la and b) our FPD 
method reasonably describe the basic behaviour of a ro¬ 
tating disk. Here we show an additional example which 
indicates the validity of our FPD method. Supplemen¬ 
tary Figure 1 shows the relation between the effective 
exerted torque Teff and the angular velocity of a disk. 
We can see a perfect linear relation between them: H = 
Teff/(47r?7ag^). Here Teff is the effective torque for our 
disk, which is described by the smooth profile function 
0(r) and given by Teff = T / |rp(/)(r)^dr/ / |rp0(r)(ir, 
where T is the exerted torque. From this analysis, we 
obtain Ueff = 1.08a. 


Supplementary Note 2: The behaviour of a pair 
of rotating discs 

As briefly discussed in the main text, a pair of rotating 
disks whose rotating axes are parallel have repulsive hy¬ 
drodynamic interactions. So the distance between them 
monotonically increases with time (see the red curve in 
Supplementary Fig. 2). However, decreasing the size 
of the third ^ dimension of the quasi-2D simulation box 
(128^ in X and y directions) from 32 to 16 leads to the 
switching of the relaxation mode of vortex from cascade 
to anti-cascade mode. Accordingly, the two rotating 
disks can have a stable rotating trajectory with a con¬ 
stant radius (see the green curve in Supplementary Fig. 
2). This is a consequence of the force balance between 
the repulsive Magnus force and the attractive hydrody¬ 
namic force. For 2D, thus, a pair of particles can form 
a stable rotating trajectory with a constant radius, as 
shown in Fig. Id. 






